#include "AnaDef.h"
#include "BaseUtils.h"
#include "TTreeStream.h"

#define EPSILON 1e-12

/*
grep static dEdxGausN.h  | grep -v return | awk '{print $3}' | sort |wc
grep dEdxGausN dEdxGausN.h  | grep -v printf | grep class -v | grep -v grep | awk -F\: '{print $3}' | sort|wc

49
50 with commented dEtadY
*/

class dEdxGausN
{
 public: 
  enum bbopt{
    kmBB = 1,
    kRawALEPH = 2,
    kLund = 16
  };

  enum anchornopt{
    kAnchorN1E4 = 1,
    kAnchorN1E5 = 2,
    kAnchorN1E6 = 4,
    kAnchorN1E9 = 8,
    kAnchorN5E3 = 16
  };

  enum TYPEID{
    kPROTON,
    kPION,
    kKAON,
    kELECTRON
  };

  //==================================================================================
  //                              controls
  //==================================================================================
 public:
  static const Char_t* Opt(){return fOpt.Data();}
  static Int_t Ntype(){CheckIni(); return fNtype;}
  static Double_t MinResolution(){CheckIni();return fMinResolution;}
  static Double_t MaxResolution(){CheckIni();return fMaxResolution;}
  static Double_t IniTolerance(){CheckIni(); return 2;}
  static Int_t IniStrategy(){return 1;}
  static Double_t IniResolution(){CheckIni(); return fIniResolution;}
  static Bool_t IsWithRef(){return fkWithRef;}
  static Bool_t IsNoHighDensity(){return fkNoHighDensity;}
  static Bool_t IsFF(){return fOpt.Contains("FF");}

  static Int_t NparType() {CheckIni(); return fNparType;}
  static BaseUtils::FFunc GausFunc(){CheckIni(); return fGausFunc;}

  static Int_t BBOpt(){return fBBopt;}
  static Int_t NparBB(){CheckIni();return fNparBB;}
  static const Double_t * IniBBpar(){CheckIni(); return fIniBBpar;}
  static BaseUtils::FFunc BBFunc(){CheckIni();return fBBFunc;}

  static Int_t LastfithIter(){return 1;}  
  static Int_t GetkETA(){return fkETA;}
  static Double_t Eta0(){return extEta0(fkETA);}
  static Double_t Eta1(){return extEta1(fkETA);}

  static Double_t extEta0(const Int_t keta){return keta*0.2;}
  static Double_t extEta1(const Int_t keta){return keta*0.2+0.2;}
 
  static Bool_t IsEtaOK(const Double_t eta);

  static Int_t GetkCEN(){return fkCEN;}
  static Bool_t IsCenOK(const Double_t cen);

  static TString GetEffFile();
  static TString GetSecoFile();
  static Int_t GetkCharge();

  static const Char_t * TOFTag(){
    const Char_t* tags[]={"","1sTOF","2sTOF"}; 
    return tags[fTOFOpt];
  }

 private:
  static Double_t BetheBlochGeant(Double_t bg,
                                  Double_t kp0=2.33,
                                  Double_t kp1=0.20,
                                  Double_t kp2=3.00,
                                  Double_t kp3=173e-9,
                                  Double_t kp4=0.49848
                                  );

  static Double_t Cen0(){if(!fCenBins){printf("null fCenBins!!\n"); exit(1);} return fCenBins[fkCEN];}
  static Double_t Cen1(){if(!fCenBins){printf("null fCenBins!!\n"); exit(1);} return fCenBins[fkCEN+1];}
  
  static Int_t LowPtSpacing(){return 1;}

  static TString fOpt;
  static Int_t fNtype;
  static Double_t fMinResolution;
  static Double_t fMaxResolution;
  static Double_t fIniResolution;

  static Int_t fNparType;
  static BaseUtils::FFunc fGausFunc;

  static Int_t fBBopt;
  static BaseUtils::FFunc fBBFunc;
  static Int_t fNparBB;
  static Double_t *fIniBBpar;

  static Int_t fAnchorOpt;
  static Bool_t fkEnergyLossCorrection;
  static Bool_t fkWithRef;
  static Bool_t fkNoHighDensity;
  static Int_t fkETA;
  static Int_t fkCEN;
  static Double_t *fCenBins;

  static Int_t fTOFOpt;

  //==================================================================================
  //                               initialization
  //==================================================================================
 public:
  static void EtaIni();
  static void CenIni();
  static void Ini();
  
 private:
  static void CheckIni();
  
  //==================================================================================
  //                                Parameters
  //==================================================================================
 public:
  static Int_t Npar(){return NparType()*Ntype();}

  static void ResetIfixed(Bool_t pfix[]);
  static void SetMean(const Double_t xx, const Double_t bbpar[], Double_t pregauspar[], Int_t *neworder=0x0, const Bool_t kprint=0);
  static void ToGausPar(const Double_t *prepar, Double_t *gauspar, const Double_t *prelowerr=0x0, const Double_t *prehigherr=0x0, Double_t *gauslowerr=0x0, Double_t *gaushigherr=0x0);

 private:  
  static void SetGausErr(Double_t * gauserr, const Double_t * preerr, const Double_t * prepar);
 
  //==================================================================================
  //                             functions
  //==================================================================================
 public: 
  static Double_t Func(const Double_t *xx, const Double_t *raw);
  static Double_t GetMean(const Double_t xx, const Double_t mass, const Double_t *bbpar, const Bool_t kprint=0);
  static Double_t Material();

 private:  
  static Double_t SatuAl(const Double_t *xx, const Double_t *par);
  static Double_t SatuBb(const Double_t *xx, const Double_t *par);
  static Double_t myBetheBloch(const Double_t *xx, const Double_t *par);
  static Double_t RawAleph(const Double_t *xx, const Double_t *par);
  static Double_t Lund(const Double_t *xx, const Double_t *par);
  static Double_t SatuLu(const Double_t *xx, const Double_t *par);
  
  //==================================================================================
  //                             algorithms
  //==================================================================================
 public:
  static void Swap(const Double_t xx, const Int_t targetip, Double_t pars[], const Int_t oldid[], const Int_t newid[]);
  static void GetInvMap(const Int_t order[], Int_t invmap[]);

 private:  
  static void GetInvMap(const Int_t nn, const Int_t order[], Int_t invmap[]);
  
  //==================================================================================
  //                               kinematics
  //==================================================================================
 public:
  static Double_t Mass(const Int_t id);
  static Int_t PDG2Type(const Int_t pdg);
  static Double_t DeltaY(const Double_t mass, const Double_t pt, const Double_t eta0, const Double_t eta1);
  static Double_t RatioEtaCorrection(const Int_t itype, const Double_t ptmean, const Double_t eta0, const Double_t eta1);
  static Double_t dYRatio(const Int_t itype, const Double_t ptmean, const Int_t etai);

  static Double_t Mproton   (){ return  0.938272;}
  static Double_t Mpion     (){ return  0.139570;}
  static Double_t Mkaon     (){ return  0.493677;}
  static Double_t Melectron (){ return  5.109989e-04;}
  
 private:  
  static Double_t Eta2Y(const Double_t mass, const Double_t pt, const Double_t eta);
  
  //==================================================================================
  //                            histogram
  //==================================================================================
 public:
  static void GetHFit(const TH2D *hrawdata, const TH1D * hrawyield, TH2D * & hfitdata, TH1D *& hfityield);
  
  //between 0.6-1 GeV p/K/pi/e are separated
  //ix=18 has too large width for proton
  //roughly corresponds to log10(0.1) ~ log10(20) bin 19
  //if fith start at 170, then 185 is well described in fithRev3
  static Int_t StartBin(){return 17* (fkNoHighDensity ? 1 : HigherDensity() )*ScaleUpNXbins();}

  //according to this, find the bin as reduction reference
  //1e4 for all
  //non gausiann tail appears at 10^-4 level
  static Double_t PbinThres(){return /*200*/ 50;}

  static Bool_t FilterIX(const Int_t ix);//, const Double_t xx);

  static TH1D *GetHistPtSpectrum();
  static TH2D* GetHistTOFVsLog10P(const TString hn);
  static TH2D *GetHistTOFMass(const TString hn);
  static TH2D *GetHistTOFpar(const Int_t ipar,const TString hn);
  static TH2D* GetHistLog10PtVsLog10P();
  static TH2D* GetHistTPCsigVsLog10P(const TString hn);
  static TH2D* GetHistDCAxyVsLog10Pt();
  static TH2D* GetHistDCAzVsLog10Pt();
  static TH2D* GetHistMCPWMVsLog10Pt(const Int_t par);
  static TH2D* GetHistMCPDGVsLog10Pt(const TString hname);
  static TH2D* GetHistMCkFailEventVsLog10Pt();
  static TH2D* GetHistMCPtVsLog10Pt();

  static TH1D* GetHyield(TH2D *h0);
  
  static Int_t Nx(){return N1x()+N2x();}
  static Double_t Xmin(){return -1; }
 
 private:  
  static Int_t ScaleUpNXbins(){return 1;}//<--- test
  static Double_t Ymin(){return 0;}
  static Double_t Ymax(){return 200;}
  static Double_t Ny(){return 400;}
 
  //X: momentum (tpcip or esd p); Y: tpcsig, pt
  //define hist, do NOT use BinLogX, becasue the bin center is not a good value for fit!!!
  //P: 0.1 ~ 100
  //log10P: -1 ~ 2
  //-1 ~ 0.6: 40 bins * 10 
  //0.6 ~ 2: 35 bins
  //xref = 0.6
  static Double_t Xmax(){return 2; }
  //< Xmid higher bin density
  static Double_t Xmid(){return 0.6;}
  static Int_t HigherDensity(){return 10;}
  
  //finer bins at low p to avoid too large bin width w.r.t. dEdx resolution (majorly protons)
  static Int_t N1x(){return 32 * HigherDensity() * ScaleUpNXbins();}//8*4 * 10
  static Int_t N2x(){return 28 * ScaleUpNXbins();}//4*7
  static void SetPPtAxis(TAxis *ax);
  
  //==================================================================================
  //                            outputs
  //==================================================================================
 public:
  static void SavePar(TTreeSRedirector * stream, const TString tn, Double_t *totalyield, Double_t *yields, Double_t *pars, Double_t *lowerrs, Double_t *higherrs, const Int_t order[], Double_t *outpar, Double_t *outlowerr, Double_t *outhigherr, Int_t *orderinv=0x0);
  static TH1D* dEdxFit(TGraphErrors * gr, const TString hname, const Bool_t kfailexit, const Double_t *inipar, const Double_t *lows=0x0, const Double_t *highs=0x0, const Double_t maxchi2=1e10);
  static void DrawTF1(const Int_t id, const TString tag,  const Double_t pars[], const Int_t col, TList *ll, FILE *log=0x0);
  static void HistPrint();

 private:  
  static void Print();

  //==================================================================================
  //                            end
  //==================================================================================

};

TString dEdxGausN::fOpt = "";
Int_t dEdxGausN::fNtype = -999;
Double_t dEdxGausN::fMinResolution = -999;
Double_t dEdxGausN::fMaxResolution = -999;
Double_t dEdxGausN::fIniResolution = -999;
Bool_t dEdxGausN::fkWithRef = kFALSE;
Bool_t dEdxGausN::fkNoHighDensity = kFALSE;

Int_t dEdxGausN::fNparType = 0;
BaseUtils::FFunc dEdxGausN::fGausFunc = 0x0;

Int_t dEdxGausN::fBBopt = -1;
BaseUtils::FFunc dEdxGausN::fBBFunc = 0x0;
Int_t dEdxGausN::fNparBB = 0;
Double_t * dEdxGausN::fIniBBpar = 0x0;

Int_t  dEdxGausN::fAnchorOpt = 0;
Bool_t dEdxGausN::fkEnergyLossCorrection = kFALSE;
Int_t dEdxGausN::fkETA=-999;
Int_t dEdxGausN::fkCEN=-999;
Double_t* dEdxGausN::fCenBins=0x0;

Int_t dEdxGausN::fTOFOpt = 0;
//==================================================================================
//                               controls
//==================================================================================
Bool_t dEdxGausN::IsEtaOK(const Double_t eta)
{
  if(fkETA==100)
    return kTRUE;

  if(TMath::Abs(eta)>Eta0() && TMath::Abs(eta)<Eta1())
    return kTRUE;

  return kFALSE;
}

Bool_t dEdxGausN::IsCenOK(const Double_t cen)
{
  if(cen>Cen0() && cen<Cen1())
    return kTRUE;

  return kFALSE;
}

//==================================================================================
//                               initialization
//==================================================================================

void dEdxGausN::EtaIni()
{
  TString pwd = gSystem->pwd();
  pwd.ToUpper();
  
  //================== fkETA
  if(pwd.Contains("ETA0"))
    fkETA = 0;
  else  if(pwd.Contains("ETA1"))
    fkETA = 1;
  else  if(pwd.Contains("ETA2"))
    fkETA = 2;
  else  if(pwd.Contains("ETA3"))
    fkETA = 3;
  else if(pwd.Contains("ETAALL"))
    fkETA=100;
  else{
    printf("dEdxGausN::EtaIni wrong eta!! %s\n", pwd.Data());exit(1);
  }
}

void dEdxGausN::CenIni()
{
  const TString pwd = gSystem->pwd();
  
  fCenBins=new Double_t[100];
  fCenBins[0]=0;
  fCenBins[1]=5;
  fCenBins[2]=10;
  fCenBins[3]=20;
  fCenBins[4]=40;
  fCenBins[5]=60;
  fCenBins[6]=80;
  fCenBins[7]=100;

  //================== fkCEN
  if(pwd.Contains("CEN0"))
    fkCEN = 0;
  else  if(pwd.Contains("CEN1"))
    fkCEN = 1;
  else  if(pwd.Contains("CEN2"))
    fkCEN = 2;
  else  if(pwd.Contains("CEN3"))
    fkCEN = 3;
  else  if(pwd.Contains("CEN4"))
    fkCEN = 4;
  else  if(pwd.Contains("CEN5"))
    fkCEN = 5;
  else  if(pwd.Contains("CEN6"))
    fkCEN = 6;
  else{
    printf("dEdxGausN::CenIni no cen!! %s\n", pwd.Data());
  }

}

Int_t dEdxGausN::GetkCharge()
{
  const TString pwd = gSystem->pwd();

  if(pwd.Contains("POS") && !pwd.Contains("NEG"))
    return 1;
  else if(pwd.Contains("NEG") && !pwd.Contains("POS"))
    return -1;
  else if(!pwd.Contains("POS") && !pwd.Contains("NEG"))
    return 0;
  else{
    printf("dEdxGausN::GetkCharge error %s\n", pwd.Data()); exit(1);
  }
}

TString dEdxGausN::GetSecoFile()
{
  Int_t kcharge = GetkCharge();
  TString dir;
  if(kcharge ==1)
    dir = "/hera/alice/xlu/workspace/hera2dbetabloch/SummarySeco_POS/";
  else if(kcharge == -1)
    dir = "/hera/alice/xlu/workspace/hera2dbetabloch/SummarySeco_NEG/";
  else
    dir = "/hera/alice/xlu/workspace/hera2dbetabloch/SummarySeco_all/";

  TString pre;
  if(fOpt.Contains("11a"))
    pre="11a";
  else if(fOpt.Contains("10d"))
    pre="10d";
  else if(fOpt.Contains("10e"))
    pre="10e";
  else if(fOpt.Contains("0de"))
    pre="10d";
  else if(fOpt.Contains("10h")){
    pre=Form("10hCEN%d", fkCEN);
  }
  else{
    printf("dEdxGausN::GetSecoFile bad fOpt %s\n",fOpt.Data());exit(1);
  }

  const TString secofile = Form("%s%s_ETA%d_seco.root",dir.Data(), pre.Data(), fkETA);
  printf("dEdxGausN::GetSecoFile %s\n", secofile.Data());
  return secofile;
}

TString dEdxGausN::GetEffFile()
{
  Int_t kcharge = GetkCharge();
  TString dir;
  if(kcharge ==1)
    dir = "/hera/alice/xlu/workspace/hera2dbetabloch/SummaryEff_POS/";
  else if(kcharge == -1)
    dir = "/hera/alice/xlu/workspace/hera2dbetabloch/SummaryEff_NEG/";
  else
    dir = "/hera/alice/xlu/workspace/hera2dbetabloch/SummaryEff_all/";

  TString pre;
  if(fOpt.Contains("11a"))
    pre="11a";
  else if(fOpt.Contains("10d"))
    pre="10d";
  else if(fOpt.Contains("10e"))
    pre="10e";
  else if(fOpt.Contains("0de"))
    pre="10d";
  else if(fOpt.Contains("10h")){
    pre=Form("10hCEN%d", fkCEN);  
  }
  else{
    printf("dEdxGausN::GetEffFile bad fOpt %s\n",fOpt.Data());exit(1);
  }

  const TString efffile = Form("%s%s_ETA%d_eff.root",dir.Data(), pre.Data(), fkETA);
  printf("dEdxGausN::GetEffFile %s\n", efffile.Data());
  return efffile;
}

void dEdxGausN::Ini()
{
  //
  //only judged by "pwd"
  //

  BaseUtils::SetMinOpt(BaseUtils::kTMINUIT);

  const TString pwd = gSystem->pwd();
  
  //==================================================================================
  //                              first level from "pwd" 
  //==================================================================================

  //================== fOpt

  if(pwd.Contains("10c"))    fOpt+=" 10c";
  if(pwd.Contains("10_de") || pwd.Contains("10-d-e"))  fOpt+=" 0de";
  if(pwd.Contains("10e"))    fOpt+=" 10e";
  if(pwd.Contains("10d"))    fOpt+=" 10d";
  if(pwd.Contains("10h"))    fOpt+=" 10h";
  if(pwd.Contains("11a"))    fOpt+=" 11a";
  if(pwd.Contains("11h"))    fOpt+=" 11h";
  if(pwd.Contains("TOY2"))   fOpt+="TOY2";
  if(pwd.Contains("TOY4"))   fOpt+="TOY4";
  if(pwd.Contains("TOY5"))   fOpt+="TOY5";

  if(fOpt.Length()>4){
    printf("dEdxGausN::Ini wrong fOpt %s %s %d\n", pwd.Data(), fOpt.Data(), fOpt.Length());exit(1);
  }
  fOpt.ReplaceAll(" ","");

  if(pwd.Contains("MC")){
    fOpt+="MC";
  }

  if(pwd.Contains("FracFunc")){
    fOpt+="FF";
  }
  //================== fBBopt
  
  fBBopt = 0;
  if(pwd.Contains("mBB"))    fBBopt += kmBB;
  if(pwd.Contains("aleph"))  fBBopt += kRawALEPH;
  if(pwd.Contains("Lund"))   fBBopt += kLund;

  //================== fAnchorOpt

  fAnchorOpt = 0;
  if(pwd.Contains("anchor3")) fAnchorOpt += kAnchorN5E3;
  if(pwd.Contains("anchor4")) fAnchorOpt += kAnchorN1E4;
  if(pwd.Contains("anchor5")) fAnchorOpt += kAnchorN1E5;
  if(pwd.Contains("anchor6")) fAnchorOpt += kAnchorN1E6;
  if(pwd.Contains("anchor9")) fAnchorOpt += kAnchorN1E9;

  //================== fTOFOpt
  fTOFOpt = 0;
  if(pwd.Contains("TOF1SIGMA")) 
    fTOFOpt = 1;
  else if(pwd.Contains("TOF2SIGMA"))
    fTOFOpt = 2;
  else
    fTOFOpt = 0;

  //================== fkEnergyLossCorrection
  //energy loss correction only derived for real data 10e, 10h, 11a!!
  if(pwd.Contains("HASENERGYLOSSCORRECTION")){
    fkEnergyLossCorrection = kTRUE;
  }
  else{
    fkEnergyLossCorrection = kFALSE;
  }

  //================== fkWithRef
  if(pwd.Contains("WithRef")){
    fkWithRef = kTRUE;
  }
  else{
    fkWithRef = kFALSE;
  }

  EtaIni();
  CenIni();

  //==================== rebin X for J pt with NoHighDensity
  if(pwd.Contains("NoHighDensity")){
    fkNoHighDensity = kTRUE;
  }
  else{
    fkNoHighDensity = kFALSE;
  }

  //==================================================================================
  //                              second level from internal control
  //==================================================================================

  //================== fNtype, MinResolution, MaxResolution, IniTolerance, IniResolution

  if(fOpt.Contains("10h")){
    fNtype = 4;
    fMinResolution = 0.03;
    fMaxResolution = 0.10;
    fIniResolution = 0.068;
  }
  else if(fOpt.Contains("10e")){
    fNtype = 4;
    fMinResolution = 0.03;
    fMaxResolution = 0.10;
    fIniResolution = 0.052;
  }
  else if(fOpt.Contains("10d")){
    fNtype = 4;
    fMinResolution = 0.03;
    fMaxResolution = 0.10;
    fIniResolution = 0.052;
  }
  else if(fOpt.Contains("0de")){
    fNtype = 4;
    fMinResolution = 0.03;
    fMaxResolution = 0.10;
    fIniResolution = 0.052;
  }
  else if(fOpt.Contains("11a")){
    fNtype = 4;
    fMinResolution = 0.03;
    fMaxResolution = 0.10;
    fIniResolution = 0.052;
  }
  else if(fOpt.Contains("10c")){
    fNtype = 4;
    fMinResolution = 0.03;
    fMaxResolution = 0.10;
    fIniResolution = 0.052;
  }
  else if(fOpt.Contains("11h")){
    fNtype = 4;
    fMinResolution = 0.03;
    fMaxResolution = 0.10;
    fIniResolution = 0.072;
  }
  else{
    printf("dEdxGausN::Ini unknown fOpt! %s\n", fOpt.Data()); exit(1);
  }

  //================== fNparType, fGausFunc

  fNparType = 3;
  fGausFunc = BaseUtils::ResolutionGaus;

  //================== fBBFunc, fNparBB, fIniBBpar

  if(fBBopt==kRawALEPH){
    fBBFunc = RawAleph;
    fNparBB = 5;
    fIniBBpar = new Double_t[fNparBB];

    const Double_t tmp[]={15, 3, log(0.03), 3, 1};
    memcpy(fIniBBpar, tmp, sizeof(Double_t)*fNparBB);   
  }
  else if(fBBopt==kmBB){
    fBBFunc = myBetheBloch;
    fNparBB = 5;
    fIniBBpar = new Double_t[fNparBB];

    const Double_t tmp[]={40, 2,  0.2, -0.2, 80};
    memcpy(fIniBBpar, tmp, sizeof(Double_t)*fNparBB);
  }
  else if(fBBopt==kLund){
    fBBFunc = Lund;
    fNparBB = 5;
    fIniBBpar = new Double_t[fNparBB];

    const Double_t tmp[]={34, 0.25, 4, 1.3, 80 };
    memcpy(fIniBBpar, tmp, sizeof(Double_t)*fNparBB);   
  }
  else{
    printf("dEdxGausN::Ini no correct BBopt specified Then set to -999!!! %d\n", fBBopt);
    fBBopt = -999;
    //printf("dEdxGausN::Ini wrong fBBopt %d\n", fBBopt);exit(1);
  }

  Print();
}

void dEdxGausN::CheckIni()
{
  if(fNtype<0 || fMinResolution <0 || fMaxResolution<0 || fIniResolution <0 || (fBBopt>0 && (!fBBFunc || fNparBB==0)) || (!fGausFunc || fNparType == 0) ){
    printf("dEdxGausN::CheckIni() fail! %d %f %f %f -- %d %d\n", fNtype, fMinResolution, fMaxResolution, fIniResolution, fNparBB, fNparType); exit(1);
  }
}

//==================================================================================
//                                Parameters
//==================================================================================

void dEdxGausN::ResetIfixed(Bool_t pfix[])
{
  for(Int_t ipar = 0; ipar<Npar(); ipar++){
    pfix[ipar] = kFALSE;
  }
}

void dEdxGausN::SetMean(const Double_t xx, const Double_t bbpar[], Double_t pregauspar[], Int_t *neworder, const Bool_t kprint)
{
  //
  //bbpar protected by "const"
  //
  //=============== set mean ==================

  //printf("longtest setmean %e %e %e %e %e\n", bbpar[0], bbpar[1], bbpar[2], bbpar[3], bbpar[4]);

  Double_t rawys[Ntype()];
  for(Int_t itype=0; itype<Ntype(); itype++){
    rawys[itype] = GetMean(xx, Mass(itype), bbpar, kprint);
    //printf("longtest %d %e\n", itype, rawys[itype]);
  }

  Int_t id[Ntype()];
  TMath::Sort(Ntype(), rawys, id, kFALSE);

  Double_t sortedys[Ntype()];
  for(Int_t ith=0; ith<Ntype(); ith++){
    sortedys[ith] = rawys[id[ith]];
    //to repulsive parametrization of [+1]
    if(ith){
      sortedys[ith] -= rawys[id[ith-1]];
    }
  }

  for(Int_t itype=0; itype<Ntype(); itype++){
    pregauspar[itype*NparType()+1] = sortedys[itype];
  }

  if(neworder){
    //update neworder
    memcpy(neworder, id, sizeof(id));
  }
}

void dEdxGausN::ToGausPar(const Double_t *prepar, Double_t *gauspar, const Double_t *prelowerr, const Double_t *prehigherr, Double_t *gauslowerr, Double_t *gaushigherr)
{
  //dummy paramters has to be set as "fixed" in TMinuit
  //no change is made to input pre*. Protected by "const"

  //============================================================
  //prepar[NparType()*itype + 0 ]: fraction
  //============================================================

  Double_t yields[Ntype()];
  for(Int_t itype=0; itype<Ntype(); itype++){
    //keep abs, this is the internal fit par used in TMinuit.
    //i.e. even fPar[0]<0, TMinuit still see a's > 0
    yields[itype] = TMath::Abs( prepar[ itype*NparType() ] );
  }

  //============================================================
  //prepar[NparType()*itype + 1 ]: pre-mean p1, mean_i = \sum_0^i-1 abs(p1)
  //============================================================

  Double_t means[Ntype()];
  Double_t mm = 0;
  for(Int_t itype = 0; itype < Ntype(); itype++){
    //all positive
    // in this way the two peaks can not overlap and the minimum would not be the averaged peak
    mm += TMath::Abs(prepar[itype*NparType()+1]);
    means[itype] = mm;
  }

  //============================================================
  //prepar[NparType()*itype + 2 ]: always resolution
  //============================================================

  Double_t widths[Ntype()];
  for(Int_t itype=0; itype<Ntype(); itype++){
    widths[itype] = prepar[ itype*NparType()+2]; 
  }

  //============================================================
  //directory copy of other parameters
  //============================================================

  for(Int_t itype = 0; itype < Ntype(); itype++){
    gauspar[itype*NparType() + 0]=yields[itype];
    gauspar[itype*NparType() + 1]=means[itype];
    gauspar[itype*NparType() + 2]=widths[itype];
  }

  if(!prelowerr || !gauslowerr)
    return;

  SetGausErr(gauslowerr, prelowerr, prepar);
  SetGausErr(gaushigherr, prehigherr, prepar);
}

void dEdxGausN::SetGausErr(Double_t * gauserr, const Double_t * preerr, const Double_t * prepar)
{
  Double_t eM2=0;//mean error ^2
  for(Int_t itype = 0; itype < Ntype(); itype++){
    eM2 += TMath::Power(preerr[1],2);
    gauserr[itype*NparType()+1] = TMath::Sqrt(eM2);

    for(Int_t ipar = 0; ipar<NparType(); ipar++){
      if(ipar==1){
        continue;
      }

      const Int_t id = itype*NparType()+ipar;
      gauserr[id] = preerr[id];
    }
  }
}

//==================================================================================
//                             functions
//==================================================================================
Double_t dEdxGausN::Material()
{
  //const Double_t material = 6.5e-3;//AliMathBase::BetheBlochAleph
  //const Double_t material = 4.06;//Solid, //Oct 16, 2012, 15:20 test, averaged among 10e, 10h, p/K
  //derived from more precise eta cut//Di 15. Jan 16:26:14 CET 2013

  //derived from 4-eta-bin scan//Mi 16. Jan 17:43:37 CET 2013
  //TF1 * fm = new TF1("fall","[0]/TMath::Sin(2*TMath::ATan(TMath::Exp(-x)))",-1,1);
  //fm->SetParameter(0,3.287829);

  //fm->Eval(0.1) //=3.30428184885474474e+00
  //fm->Eval(0.3) //=3.43689428156977694e+00
  //fm->Eval(0.5) //=3.70744134955852944e+00
  //fm->Eval(0.7) //=4.12678105661457728e+00

  Double_t material = -999;
  if(fkETA==0)//eta 0.0-0.2
    material = 3.30428;
  else if(fkETA==1)
    material = 3.43689;
  else if(fkETA==2)
    material = 3.70744;
  else if(fkETA==3)//eta 0.6-0.8
    material = 4.12678;
  else{
    printf("dEdxGausN::Material wrong eta and material!! %d %f\n", fkETA, material);exit(1);
  }
 
  return material; 
}

//hard coded material!!
//http://alisoft.cern.ch/viewvc/trunk/STEER/STEERBase/AliExternalTrackParam.cxx?view=markup&root=AliRoot
Double_t dEdxGausN::BetheBlochGeant(Double_t bg,
                                                Double_t kp0,
                                                Double_t kp1,
                                                Double_t kp2,
                                                Double_t kp3,
                                                Double_t kp4) {
  //
  // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
  //
  // bg  - beta*gamma
  // kp0 - density [g/cm^3]
  // kp1 - density effect first junction point
  // kp2 - density effect second junction point
  // kp3 - mean excitation energy [GeV]
  // kp4 - mean Z/A
  //
  // The default values for the kp* parameters are for silicon. 
  // The returned value is in [GeV/(g/cm^2)].
  // 

  const Double_t mK  = 0.307075e-3; // [GeV*cm^2/g]
  const Double_t me  = 0.511e-3;    // [GeV/c^2]
  const Double_t rho = kp0;
  const Double_t x0  = kp1*2.303;
  const Double_t x1  = kp2*2.303;
  const Double_t mI  = kp3;
  const Double_t mZA = kp4;
  const Double_t bg2 = bg*bg;
  const Double_t maxT= 2*me*bg2;    // neglecting the electron mass
  
  //*** Density effect
  Double_t d2=0.; 
  const Double_t x=TMath::Log(bg);
  const Double_t lhwI=TMath::Log(28.816*1e-9*TMath::Sqrt(rho*mZA)/mI);
  if (x > x1) {
    d2 = lhwI + x - 0.5;
  } else if (x > x0) {
    const Double_t r=(x1-x)/(x1-x0);
    d2 = lhwI + x - 0.5 + (0.5 - lhwI - x0)*r*r*r;
  }

  return mK*mZA*(1+bg2)/bg2*
    (0.5*TMath::Log(2*me*bg2*maxT/(mI*mI)) - bg2/(1+bg2) - d2);
}


Double_t dEdxGausN::GetMean(const Double_t xx, const Double_t mass,  const Double_t *bbpar, const Bool_t kprint)
{
  Double_t pp = TMath::Power(10,xx); //xx=log10(p)

  if(fkEnergyLossCorrection){
    const Double_t material = Material();
    pp = TMath::Sqrt(TMath::Power(TMath::Sqrt(pp*pp+mass*mass)-material*BetheBlochGeant(pp/mass),2)-TMath::Power(mass,2));
  }

  const Double_t bg = pp/mass;

  if(kprint){
    printf("getmean %e %e %e --- %e %e %e %e %e\n", TMath::Power(10,xx), pp, bg, bbpar[0], bbpar[1], bbpar[2], bbpar[3], bbpar[4]);
  }
  //with +EPSILON to protect mean = 0 -> resolution -> inf
  //Max() to reject negative mean
  return TMath::Max(EPSILON, BBFunc()(&bg, bbpar));
}

Double_t dEdxGausN::Func(const Double_t *xx, const Double_t *inraw)
{
  Double_t par[Npar()];
  ToGausPar(inraw, par);

  return BaseUtils::GausN(xx, par, Ntype(), GausFunc(), NparType());
}

Double_t dEdxGausN::SatuBb(const Double_t *xx, const Double_t *par)
{
  const Double_t sig = TMath::Max(EPSILON, myBetheBloch(xx, par));
  return sig*BaseUtils::ProtectedExp(par[5]/sig);
}

Double_t dEdxGausN::myBetheBloch(const Double_t *xx, const Double_t *par)
{
  //RawALEPH has numerical difficulty in ln(p2+1/pow(bg, p4))
  //use flat dedx for electron. deviation from this at high bg is not inportant because anyway there the statistics is very limited
  //npar = 5;

  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = par[0];//overall prefactor
  const Double_t p1 = par[1];//beta exponent ~ 2
  const Double_t p2 = par[2];
  const Double_t p3 = par[3];
  const Double_t p4 = par[4];//plateau ~ 80

  //EPSILON added to protect 0 denominator!
  const Double_t bethebloch = p0 *( 1 + p2*TMath::Log(bg) + p3*TMath::Log(beta) )/( BaseUtils::ProtectedPow(beta, p1)+EPSILON) ;

  //tp = 100, electron momentum = 51 MeV, pion momentum = 14 GeV
  const Double_t tp = 10;

  Double_t qq = 0;
  if(bg<tp){
    qq = bethebloch;
  }
  else{
    qq = TMath::Min(p4, bethebloch);
  }

  return qq;
}

Double_t dEdxGausN::SatuAl(const Double_t *xx, const Double_t *par)
{
 
  const Double_t sig = TMath::Max(EPSILON, RawAleph(xx, par));
  return sig*BaseUtils::ProtectedExp(par[5]/sig);
}

Double_t dEdxGausN::RawAleph(const Double_t *xx, const Double_t *par)
{
  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = TMath::Abs(par[0]);
  const Double_t p1 = TMath::Abs(par[1]);
  const Double_t p2 = par[2];//TMath::Abs(par[2]);//still first have to be positive
  const Double_t p3 = TMath::Abs(par[3]);
  const Double_t p4 = TMath::Abs(par[4]);//has to be positive so that in bg-> inf there is no bg dependence anymore

  //return p0*((p1-TMath::Log(p2+1/BaseUtils::ProtectedPow(bg,p4)))/BaseUtils::ProtectedPow(beta,p3) - 1 );
  //return p0*( (p1- TMath::Log(p2+EPSILON)- TMath::Log(1+1/(BaseUtils::ProtectedPow(bg,p4)*p2+EPSILON))  )/( BaseUtils::ProtectedPow(beta,p3)+EPSILON ) - 1 );
  //log(p2) -> p2 because p2 sometimes ~EPSILON
  return p0*( (p1- p2 - TMath::Log(1+1/(BaseUtils::ProtectedPow(bg,p4)*BaseUtils::ProtectedExp(p2)+EPSILON))  )/( BaseUtils::ProtectedPow(beta,p3)+EPSILON ) - 1 );
  //p2 can be very small ~ 1e-9, do not separate p2 out! -> absorb p2 into p1, remember to change initial value of p1
  //return p0*( (p1- TMath::Log(1+1/(BaseUtils::ProtectedPow(bg,p4)*p2+EPSILON))  )/( BaseUtils::ProtectedPow(beta,p3)+EPSILON ) - 1 );
  //change p0 as plateau
  //return p0/(p1-1)*( (p1- TMath::Log(1+1/(BaseUtils::ProtectedPow(bg,p4)*p2+EPSILON))  )/( BaseUtils::ProtectedPow(beta,p3)+EPSILON ) - 1 );
  //redefine p1 so that p1>1 --> p1>0
  //return p0/p1*( (p1 + 1 - TMath::Log(1+1/(BaseUtils::ProtectedPow(bg,p4)*p2+EPSILON))  )/( BaseUtils::ProtectedPow(beta,p3)+EPSILON ) - 1 );
}

Double_t dEdxGausN::SatuLu(const Double_t *xx, const Double_t *par)
{
  const Double_t sig = TMath::Max(EPSILON, Lund(xx, par));
  return sig*BaseUtils::ProtectedExp(par[5]/sig);
}

Double_t dEdxGausN::Lund(const Double_t * xx,  const Double_t * rawpar)
{
  Double_t par[100];
  memcpy(&(par[1]), rawpar, sizeof(Double_t)*5);

  const Double_t bg = xx[0];

  // bg is beta*gamma
  const Double_t beta2 = bg*bg / (1.0+bg*bg);
  const Double_t a = par[1];
  const Double_t b = par[2]*par[1];//in this way value = rawpar[0]*()
  const Double_t c = par[3];
  const Double_t e = par[4];
  const Double_t f = par[5];
  const Double_t d = TMath::Exp(c*(a-f)/b);

  const Double_t powbg = TMath::Power(1.0+bg, c);

  const Double_t value = a/TMath::Power(beta2,e) +
    b/c*TMath::Log(powbg/(1.0 + d*powbg));
  return value;
}

//==================================================================================
//                             algorithms
//==================================================================================
void dEdxGausN::GetInvMap(const Int_t order[], Int_t invmap[])
{
  GetInvMap(Ntype(), order, invmap);
}

void dEdxGausN::GetInvMap(const Int_t nn, const Int_t order[], Int_t invmap[])
{
 for(Int_t itype=0; itype<nn; itype++){
    invmap[itype]=-999;
    for(Int_t ss=0; ss<nn; ss++){
      if(order[ss]==itype){
        invmap[itype]=ss;
        break;
      }
    }
    if( invmap[itype]<0 ){
      printf("dEdxGausN::GetInvMap fail to find invmap %d!\n", itype);exit(1);
    }
  }
}
void dEdxGausN::Swap(const Double_t xx, const Int_t targetip, Double_t pars[], const Int_t oldid[], const Int_t newid[])
{
  Double_t oldpar[Npar()];
  memcpy(oldpar, pars, sizeof(oldpar));
  
  for(Int_t newith=0; newith<Ntype(); newith++){
    Int_t oldith = -999;
    for(Int_t ss=0; ss<Ntype(); ss++){
      if(oldid[ss]==newid[newith]){  //the newith-th par now = the ss-th par before
        oldith=ss;
        break;
      }
    }
    if(oldith<0){
      printf("dEdxGausN::Swap old id not found!! error!!\n");exit(1);
    }
    
    if(newith!=oldith){
      printf("dEdxGausN::Swap reordering %d at xx %e: type %d old %d new %d\n", targetip, xx, newid[newith], oldith, newith);
     
      pars[NparType()* newith + targetip] = oldpar[NparType()* oldith +targetip];
    }
  }
  
}

//==================================================================================
//                               kinematics
//==================================================================================

Double_t dEdxGausN::dYRatio(const Int_t itype, const Double_t ptmean, const Int_t etai)
{
  //
  //dy_i / dy_0, 0: mid rapidity
  //
  const Double_t y0 = DeltaY(Mass(itype), ptmean, extEta0(0), extEta1(0));
  const Double_t yi = DeltaY(Mass(itype), ptmean, extEta0(etai), extEta1(etai));
  return yi/y0;
}

Double_t dEdxGausN::RatioEtaCorrection(const Int_t itype, const Double_t ptmean, const Double_t eta0, const Double_t eta1)
{
  //
  //eta correction for yield ratio of Nitype/Npion, should correct N -> N/dy
  //

  const Double_t dypion = DeltaY(Mass(kPION), ptmean, eta0, eta1);
  const Double_t dytype = DeltaY(Mass(itype), ptmean, eta0, eta1);
  if(dytype<EPSILON){
    printf("dEdxGausN::RatioEtaCorrection dytype<eps ! %d %e %e %e %e\n", itype, ptmean, eta0, eta1, dytype);exit(1);
  }
  return dypion/dytype;
}

Double_t dEdxGausN::Mass(const Int_t id)
{
  switch(id){
  case kPROTON:
    return Mproton();
  case kPION:
    return Mpion();
  case kKAON:
    return Mkaon();
  case kELECTRON:
    return Melectron();
  default:
    printf("dEdxGausN::Mass wrong id %d\n", id); exit(1);
  }

  return -999;
}

Int_t dEdxGausN::PDG2Type(const Int_t pdg)
{
  Int_t itype = -999;
  switch(pdg){
  case 11:
    itype = kELECTRON;
    break;
  case 211:
    itype = kPION;
    break;
  case 2212:
    itype = kPROTON;
    break;
  case 321:
    itype = kKAON;
    break;
  default:
    itype = -999;
    break;
  }
  return itype;
}

Double_t dEdxGausN::Eta2Y(const Double_t mass, const Double_t pt, const Double_t eta)
{
  //checked!
  /*
    TLorentzVector b
    b.SetPtEtaPhiM(1, 0.7, 0, 1)
    double eta = 0.7, mass = 1, pt = 1
    const Double_t sinh = TMath::SinH(eta)/TMath::Sqrt(1+TMath::Power(mass/pt,2));
    TMath::ASinH(sinh) - b.Rapidity()
    b.SetPtEtaPhiM(1, 0.7, 3, 1)
    TMath::ASinH(sinh) - b.Rapidity()
   */

  const Double_t sinhy = TMath::SinH(eta)/TMath::Sqrt(1+TMath::Power(mass/pt,2));
  return TMath::ASinH(sinhy);
}

Double_t dEdxGausN::DeltaY(const Double_t mass, const Double_t pt, const Double_t eta0, const Double_t eta1)
{
  const Double_t y0 = Eta2Y(mass, pt, eta0);
  const Double_t y1 = Eta2Y(mass, pt, eta1);
  return TMath::Abs(y0-y1);
}

/*
Double_t dEdxGausN::dEtady(const Double_t mass, const Double_t pp)
{
  //differential
  const Double_t ee = TMath::Sqrt(mass*mass+pp*pp);
  return ee/pp;
}
*/

//==================================================================================
//                            histogram
//==================================================================================

void dEdxGausN::SetPPtAxis(TAxis *ax)
{
  const Double_t dx1 = (Xmid()-Xmin())/N1x();
  const Double_t dx2 = (Xmax()-Xmid())/N2x();

  Double_t bins[Nx()+1];
  bins[0] = Xmin();
  //[Xmin, Xmid] N1x bins
  for(Int_t ibin=1; ibin<=N1x(); ibin++){
    bins[ibin] = bins[ibin-1]+dx1;
  }
  //[Xmid, Xmax] N2x bins
  for(Int_t ibin=1; ibin<=N2x(); ibin++){
    bins[N1x()+ibin] = bins[N1x()+ibin-1]+dx2;
  }

  ax->Set(Nx(), bins);

  printf("\n");
  for(Int_t ibin=0; ibin<=ax->GetNbins(); ibin++){
    const Double_t xh = ax->GetBinUpEdge(ibin);
    printf("dEdxGausN::SetPPtAxis %d: %e = log10 %e\n", ibin, xh, TMath::Power(10, xh));
  }
  printf("\n");
}

TH1D *dEdxGausN::GetHistPtSpectrum()
{
  //bin width has to be small in low pt, otherwise there is delta pt smearing
  //SetPPtAxis (log10 -1 ~ 2 in 60 bins)
  const Int_t nbin = 600;
  Double_t bins[nbin+1]={0.1};
  const Double_t dlogx = 0.05/10.;

  for(Int_t ii=1; ii<=nbin; ii++){
    bins[ii]=bins[ii-1]*TMath::Power(10,dlogx);
  }

  TH1D * hh = new TH1D("ptSpectrum","pt", nbin, bins[0], bins[nbin]);

  hh->GetXaxis()->Set(nbin, bins);


  //no!!!! pt spectrum is everywhere steeply falling
  /*
  //not too small, otherwise too large stat. error
  TH1D * hh = new TH1D("ptSpectrum","pt",  Nx(), Xmin(), Xmax());
  SetPPtAxis(hh->GetXaxis());
  */
  return hh;
}

TH2D* dEdxGausN::GetHistTOFVsLog10P(const TString hn)
{
  TH2D * hh = new TH2D(hn,"tof id : log10 p", Nx(), Xmin(), Xmax(), 4, 0, 4);
  SetPPtAxis(hh->GetXaxis());
  return hh;
}

TH2D *dEdxGausN::GetHistTOFMass(const TString hn)
{
  TH2D * hh = new TH2D(hn,"tof mass : log10 p", Nx(), Xmin(), Xmax(), 200, 0,2);
  SetPPtAxis(hh->GetXaxis());
  return hh;
}

TH2D *dEdxGausN::GetHistTOFpar(const Int_t ipar, const TString hn)
{
  TH2D * hh = new TH2D(Form("%s%d",hn.Data(), ipar), Form("tof %d nsigmas : log10 p", ipar), Nx(), Xmin(), Xmax(), 400, -20, 20);
  SetPPtAxis(hh->GetXaxis());
  return hh;
}

TH2D* dEdxGausN::GetHistMCPtVsLog10Pt()
{
  TH2D * hh = new TH2D("mcpt", "mc pt : log10 p (or mom)", Nx(), Xmin(), Xmax(), Nx(), Xmin(), Xmax());
  SetPPtAxis(hh->GetXaxis());
  return hh;
}

TH2D* dEdxGausN::GetHistMCPWMVsLog10Pt(const Int_t par)
{
  TH2D * hh = new TH2D(par<4? Form("mcpwm%d",par) : "mcpwm" , Form("mc primary/weak/material par %d : log10 pt", par), Nx(), Xmin(), Xmax(), 8, 0, 8);
  SetPPtAxis(hh->GetXaxis());
  return hh;
}

TH2D* dEdxGausN::GetHistDCAxyVsLog10Pt()
{
  TH2D * hh = new TH2D("dcaxyVsPt","dcaxy : log10 pt", Nx(), Xmin(), Xmax(), 200, -0.7, 0.7);
  SetPPtAxis(hh->GetXaxis());
  return hh;
}

TH2D* dEdxGausN::GetHistDCAzVsLog10Pt()
{
  TH2D * hh = new TH2D("dcazVsPt","dcaz : log10 pt", Nx(), Xmin(), Xmax(), 200, -3, 3);
  SetPPtAxis(hh->GetXaxis());
  return hh;
}

TH2D* dEdxGausN::GetHistMCPDGVsLog10Pt(const TString hname)
{
  TH2D * hh = new TH2D(hname,"mcpdg : log10 pt", Nx(), Xmin(), Xmax(), 4, 0, 4);
  SetPPtAxis(hh->GetXaxis());
  return hh;
}

TH2D * dEdxGausN::GetHistMCkFailEventVsLog10Pt()
{
  TH2D * hh = new TH2D("kFailEventVsPt", "kFailEvent : log10 pt", Nx(), Xmin(), Xmax(), 5, 0, 5);
  SetPPtAxis(hh->GetXaxis());
  return hh;
}

TH2D* dEdxGausN::GetHistLog10PtVsLog10P()
{
  //very fine pt bins to ensure accurate mean and rms calculation
  TH2D * kk = new TH2D("kk","log10 pt : log10p (inequal x bin width)", Nx(), Xmin(), Xmax(), 10000, Xmin(), Xmax());

  SetPPtAxis(kk->GetXaxis());

  return kk;
}

TH2D* dEdxGausN::GetHistTPCsigVsLog10P(const TString hn)
{
  TH2D * hh = new TH2D(hn, hn+ " tpcsig : log10p (inequal x bin width)", Nx(), Xmin(), Xmax(), Ny(), Ymin(), Ymax());

  SetPPtAxis(hh->GetXaxis());

  return hh;
}

Bool_t dEdxGausN::FilterIX(const Int_t ix)
{

  //return 1 if remove, 0 if keep

  //to avoid the region where proton is far away which aleph is not pefect
  //if fith start at 170, then 185 is well described in fithRev3
  const Int_t newstart = StartBin();//+15;

  if(ix<newstart)
    return kTRUE;

  if(fkNoHighDensity)
    return kFALSE;

  if( ix<N1x() ){
    if((ix% (LowPtSpacing()*HigherDensity()*ScaleUpNXbins()) ) != (newstart% (LowPtSpacing()*HigherDensity()*ScaleUpNXbins()) )){
      return kTRUE;
    }
  }
  else{
    if((ix% ScaleUpNXbins() ) != (newstart% ScaleUpNXbins() )){
      return kTRUE;
    }
  }
  
  return kFALSE;
}

TH1D * dEdxGausN::GetHyield(TH2D *h0)
{
  //only consider the dEdx in range
  return h0->ProjectionX(Form("%syield", h0->GetName()), 1, h0->GetNbinsY(),"oe");
}

void dEdxGausN::GetHFit(const TH2D *hrawdata, const TH1D * hrawyield, TH2D * & hfitdata, TH1D *& hfityield)
{
  Double_t anchorN = 1e9;
  if(!fOpt.Contains("TOY")){
    if(fAnchorOpt==kAnchorN5E3)   
      anchorN = 5e3;
    else if(fAnchorOpt==kAnchorN1E4)   
      anchorN = 1e4;
    else if(fAnchorOpt==kAnchorN1E5)   
      anchorN = 1e5;
    else if(fAnchorOpt==kAnchorN1E6)   
      anchorN = 1e6;
    else if(fAnchorOpt==kAnchorN1E9)   
      anchorN = 1e9;
    else {
      printf("dEdxGausN::GetHFit wrong anchor opt!! %d\n", fAnchorOpt);exit(1);
    }
  }

  //if anchor too small, use the least N of all bins
  Double_t minn = 1e10;
  for(Int_t ix = 1; ix<=hrawyield->GetNbinsX(); ix++){
    if(FilterIX(ix)){
      continue;
    }
    
    const Double_t cc = hrawyield->GetBinContent(ix);
    if(cc<PbinThres()){
      continue;
    }
    if(cc<minn)
      minn = cc;
  }
  if(anchorN<minn){
    anchorN = minn;
    printf("check dEdxGausN::GetHFit --> anchorN modified as minimum stat %e\n", anchorN);
  }

  hfityield = (TH1D*) hrawyield->Clone("hfityield");
  hfityield->Scale(0);
  
  //rebined hrawdata!! not save in list
  TH2D *htmpori = (TH2D*) hrawdata->Clone("htmpori");
  //just include dedx in range for ToPDF
  htmpori->GetYaxis()->SetRange(1,htmpori->GetNbinsY());
  //change Ny=400 to Ny=200
  //y=100, dy=1 => dy/y = 1e-2, not enough
  //---> effect not obvious
  htmpori->RebinY(2);

  hfitdata = BaseUtils::NormalHist(htmpori, PbinThres());
  hfitdata->SetName("hfitdata");
  delete htmpori;
  
  //do not use hreduce, but just inflate horiginal to preserve more information
  //inflate hfitdata error to Anchor
  for(Int_t ix=1; ix<=hfitdata->GetNbinsX(); ix++){
    const Double_t tmpnor = hrawyield->GetBinContent(ix);
    if(tmpnor>anchorN){
      const Double_t fscale = anchorN/tmpnor;
      
      for(Int_t iy=1; iy<=hfitdata->GetNbinsY(); iy++){
        const Double_t tmperr = hfitdata->GetBinError(ix,iy);
        //sqrt important
        //so the norm in BinnedLikelihood fit is:
        //pow(A/ (sqrt(A/N)/sqrt(C/N)) ,2) /A = C
        hfitdata->SetBinError(ix,iy,tmperr/TMath::Sqrt(fscale));
      }
      
      hfityield->SetBinContent(ix, anchorN);
    }
    else{
      hfityield->SetBinContent(ix, tmpnor);
    }
  }
}

//==================================================================================
//                            outputs
//==================================================================================
void dEdxGausN::HistPrint()
{
  printf("Xmid(): %e = log10 %e\n", Xmid(), TMath::Power(10, Xmid()));

  printf("Xmin(): %e\n", Xmin());
  printf("Xmax(): %e\n", Xmax());
  printf("N1x(): %d\n", N1x());
  printf("N2x(): %d\n", N2x());
  printf("Nx(): %d\n", Nx());
  printf("Ymin(): %e\n", Ymin());
  printf("Ymax(): %e\n", Ymax());
  printf("Ny(): %e\n", Ny());
  printf("HigherDensity(): %d\n", HigherDensity());
  printf("ScaleUpNXbins(): %d\n", ScaleUpNXbins());
  printf("fkETA %d eta0 %.1f eta1 %.1f\n", fkETA, Eta0(), Eta1());
}

void dEdxGausN::Print()
{
  printf("\ndEdxGausN::Print()\n");

  printf("important controls:\n");

  printf("fBBopt %d\n", fBBopt);
  printf("fAnchorOpt: %d\n", fAnchorOpt);
  printf("fTOFOpt: %d\n", fTOFOpt);
  printf("IsWithRef(): %d;\n", IsWithRef());
  printf("IsNoHighDensity(): %d;\n", IsNoHighDensity());
  printf("fkEnergyLossCorrection: %d\n", fkEnergyLossCorrection);
  if(fkETA!=100){
    printf("Material(): %f;\n", Material());
  }
  printf("\n");

  printf("fOpt=%s; \n", fOpt.Data());
  printf("Ntype(): %d\n", Ntype());
  printf("MinResolution(): %e\n", MinResolution());
  printf("MaxResolution(): %e\n", MaxResolution());
  if(MinResolution()<0 || MaxResolution()<0){
    printf("dEdxGausN::Print min/max res < 0 !\n"); exit(1);
  }

  printf("IniResolution(): %e\n", IniResolution());

  printf("IniTolerance(): %e\n", IniTolerance());
  printf("IniStrategy(): %d\n", IniStrategy());

  printf("NparType(): %d\n", NparType());

  printf("NparBB(): %d\n", NparBB());

  printf("StartBin(): %d\n", StartBin());

  printf("PbinThres(): %e\n", PbinThres());

  printf("LowPtSpacing(): %d\n", LowPtSpacing());

  HistPrint();
  printf("\n");
}

TH1D *dEdxGausN::dEdxFit(TGraphErrors * gr, const TString hname, const Bool_t kfailexit, const Double_t *inipar, const Double_t *lows, const Double_t *highs, const Double_t maxchi2)
{
  printf("dEdxGausN::dEdxFit %s start ...\n", hname.Data());

  Double_t par[100];
  if(!inipar){
    inipar = fIniBBpar;
  }
  memcpy(par, inipar, sizeof(Double_t)* NparBB());
  Double_t err[100], chi[100];

  const Int_t kfail = BaseUtils::ChisquareFit(gr->GetN(), gr->GetX(), gr->GetY(), 0x0, gr->GetEY(), BBFunc(), NparBB(), par,err,chi, 0x0, 0x0, lows, highs); 
 
  if(kfail || chi[0]/chi[1] > maxchi2){
    printf("dEdxGausN::dEdxFit kfail %d %s or bad chi2 %e %e, %e %e\n\n", kfail,kfailexit?"exit!!!":"", chi[0], chi[1], chi[0]/chi[1], maxchi2);
    if(kfailexit){
      printf("exit!!!\n");
      exit(1);
    }
  }

  if(fBBopt==kRawALEPH){
    if(par[0]>100){
      printf("dEdxGausN::dEdxFit special check for p0 aleph, larger than 100, fall back to inipar! %e -- %e\n", par[0], inipar[0]);
      memcpy(par, inipar, sizeof(Double_t)* NparBB());
    }
  }

  TH1D *hfit = BaseUtils::GetHfit(hname, BBFunc(), par, 0.1, 4e4, kTRUE);
  TString htit;
  for(Int_t ii=0; ii<NparBB(); ii++){
    htit+=Form("%e , ",par[ii]);
  }
  hfit->SetTitle(htit);
  return hfit;
}

void dEdxGausN::DrawTF1(const Int_t id, const TString tag,  const Double_t pars[], const Int_t col, TList *ll, FILE *log)
{
  Char_t buff[2000];
  sprintf(buff, "dEdxGasuN::DrawTF1 longtest %s %d %e %e %e\n", tag.Data(), id, pars[0], pars[1], pars[2]);
  if(log){
    fprintf(log, buff,"");
  }
  else{
    printf(buff,"");
  }

  TH1D * ff = BaseUtils::GetHfit(Form("%sf%d", tag.Data(), id), GausFunc(), pars, Ymin(), Ymax());
  ff->SetLineColor(col);
  //ff->SetDrawOption("same hist C");
  ll->Add(ff);
}

void dEdxGausN::SavePar(TTreeSRedirector * stream, const TString tn,  Double_t *totalyield,  Double_t *yields, Double_t *pars,  Double_t *lowerrs,  Double_t *higherrs, const Int_t order[], Double_t *outpar, Double_t *outlowerr, Double_t *outhigherr, Int_t *orderinv)
{
  //
  //totalyield, pars, lowerrs, higherrs should not be changed
  //due to streamer, no protection is made using "const"
  //

  (*stream)<<tn<<
    "totn="<<*totalyield;

  //========================================================================
  //========================================================================

  for(Int_t itype=0; itype<Ntype(); itype++){
    for(Int_t ipar = 0; ipar<NparType(); ipar++){
      (*stream)<<tn<<
        Form("rawp%d%d=",itype, ipar)<<pars[itype*NparType()+ipar]<<
        Form("rawle%d%d=",itype, ipar)<<lowerrs[itype*NparType()+ipar]<<
        Form("rawhe%d%d=",itype, ipar)<<higherrs[itype*NparType()+ipar];
    }
  }

  //========================================================================
  //                       the following only change outpar
  //========================================================================

  //outpar and outerr have to be define outside so to have a fixed memory address because the streamer does not end here. It uses the memory address after this function.
  //pars needs to be used for next fit, can not be touched!

  ToGausPar(pars, outpar, lowerrs, higherrs, outlowerr, outhigherr);

  //========================================================
  //                         save
  //========================================================

  //the inverse map of order: invmap[itype] = ith; order[ith]=itype
  Int_t invmap[Ntype()];
  GetInvMap(order, invmap);

  for(Int_t itype=0; itype<Ntype(); itype++){
    Double_t &iyield = yields[invmap[itype]];
    iyield = *totalyield*outpar[invmap[itype]*NparType()];

    (*stream)<<tn<<
      Form("y%d=",itype)<<iyield;

    for(Int_t ipar = 0; ipar<NparType(); ipar++){
      //-> the oder of streaming have to stay always the same
      //-> save in order of particle type

      //only address of outpar and outerr can be used for streamer!
      Double_t &var = outpar[invmap[itype]*NparType()+ipar];
      Double_t &lowerr = outlowerr[invmap[itype]*NparType()+ipar];
      Double_t &higherr= outhigherr[invmap[itype]*NparType()+ipar];
      
      (*stream)<<tn<<
        Form("p%d%d=",itype, ipar)<<var<<
        Form("le%d%d=",itype, ipar)<<lowerr<<
        Form("he%d%d=",itype, ipar)<<higherr;
    }
  }

  if(orderinv){
    memcpy(orderinv, invmap, sizeof(invmap));
  }
}
